^D-AQ89  937  BRONX  COMMUNITY  COLL  NY  OEPT  OF  PHYSICS  F/6  20/3 

TRANSMISSION  OF  ELECTROMASNETIC  RADIATION  THR0U6H  NON-PLANAR  LA— ETC(U> 
MAR  80  L  A  OEACETIS  AFOSR-79-0035 

AFOSR-TR-80-0757  NL 


UNCLASSIFIED 


MICROCOPY  resolution  test  chart 

NAIIONAl  BURtAU  Of  SIANDAROS  A 


AF0Sft-1<V-'6C>35 

jjjjj/TC-80 "->7  57  I 


TRANSMISSION  OF  ELECTROMAGNETIC  RADIATION 
THROUGH  NON- PLANAR  LAYERED  MEDIA  e 


PlfiZ/fL 


1 1  !  >  ’  ■  ■'  .  7 
(£?)/  . 


/  Louis  A./DeAcetis 


<  / 


Department  of  Physics 
Bronx  Community  College-  -- 
of  The  City  University  of  New  York 
Bronx,  New  York  10453 


March- 4980  / 


Prepared  for 

Air  Force  Office  of  Scientific  Research 
Bolling  AFB 
D.C.  20332 


<?///-<=■ 

js-sssKr-' 
80“}  2  0  63 


_ liiNCUSSJFJF.Q.  _ 

SECURITY  CL ASSiFiC OF  This  PAGE  (**• n  Dmtm  tmrnrad) 


REPORT  DOCUMENTATION  PAGE 

READ  INSTRUCTIONS 

BEFORE  COMPLETING  FORM 

}  REPORT  NjMStR  ^  2.  GOVT  ACCESSION  NO. 

AFOSR.TR.  8  0-0757  /)D-A 569  9  3 

1  RECIPIENT'S  CATALOG  NUMBER 

7 

4.  TITLE  ('•rid  Subilrle) 

Transmission  of  Electromagnetic  Radiation 
through  Non-Planar  Layered  Media  / 

S  TYPE  OF  REPORT  k  PERIOD  COVERED 

Final 

«  PERFORMING  ORG.  REPORT  NUMBER 

7  AiJ  T hOR'IJ 

Louis  A.  DeAcetis 

8  CONTRACT  OR  GRANT  NUMBER',) 

AF0SR  79-0035 

9  PERFORMING  OBSAni  l  A  TiOn  NAME  and  address 

Physics  Department 

Bronx  Community  College 

Bronx,  New  York  10453 

io.  ppogp am  element. project,  task 

AREA  t  *OR<  UNIT  NUMBERS 

61 102F  2301 /A6 

II.  CONTROLLING  OFFICE  NAME  AND  ADOREsj 

AF0SR/NP 

Bolling  AFB,  Bldg.  4l0 

Washington,  DC  20332  . 

11  report  date 

March  1980  ^ 

t«  mOniTO°.nG  AGEnCy  NAME  A  ADDRESS  (it  dUlttmnt  from  Controlling  Oltico) 

IS  SECURITY  class  fol  1  hi.  r.por  1) 

UNCLASSIFIED 

i«a  DECLASSIFICATION  DOWNGRADING 
SCHEDULE 

OiSTRiBltiON  STATEMENT  (of  thf»  Rmpott) 

Approved  for  publio  release ; 
distribution  unlimited. 

yy  D»STHi0uTiOn  STATEMENT  (of  th»  •bttrmct  enrered  In  Block  30.  It  diftoront  from  Report) 


’•  Supplement  ary  notes 


1<*  KEY  WOPOS  'Conf/no#  on  t0vr$»  a»de  If  n#c#»««ry  And  identity  by  block  nutrbor) 

BESSEL  FUNCTIONS  LAYERED  MEDIA 

COMPUTER  SUBPROGRAMS  MATHEMATICAL  FUNCTIONS 

ELECTROMAGNETIC  WAltfS  RA DOMES 


*0  ABSTRACT  'Continue  aft  rotorsa  a>d#  1/  Memory  on d  tdoniHy  by  block  number,) 

An  investigation  of  transmission  of  electromagnetic  radiation 
through  non-planar  media  consisting  of  several  layers  of  dielectric 
has  begun.  A  library  search  of  recent  theoretical  and  experimental 
work  in  this  area  was  conducted,  and  a  theoretical  study  of  the 
exact  solutions  to  transmission  through  layered  cylindrical  and 
spherical  media  is  underway.  Computer  programs  to  generate  the 
spherical  and  half-integer  Bessel  Functions  of  both  the  first  and 
I  ao^r.H  1H  nd. are -Included  -in  thi.3  gSSQL^  .  ■■  ■  ..  .  - - 


DD  1473  toi’iON  or  «  hov  **  u  obsolete 


IVUSCtFlFD 


ABSTRACT 


An  investigation  of  transmission  of  electromagnetic 

radiation  through  non-planar  media  consisting  of  several 

layers  of  dielectric  has  begun.  A  library  search  of  recent 
theoretical  and  experimental  work  in  this  area  was  conducted, 

and  a  theoretical  study  of  the  exact  solutions  to  transmission 
through  layered  cylindrical  and  spherical  media  is  underway. 

The  library  survey  did  not  uncover  any  systemic*experimental 
investigations  in  this  area.  Attempts  at  expanding  exact 
solutions  to  obtain  suitable  approximations  to  the  plane- 
layer  cases  have  thus  far  not  been  successful.  Several  com¬ 
puter  programs  not  generally  available  were  developed  to 
generate  the  spherical  and  half-integer  Bessel  Functions  of 
both  the  first  and  second  kind  and  are  included  in  this  report. 
Calculations  using  these  functions  are  currently  underway. 
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introduction 


The  investigation  of  transmission  of  electromagnetic  radia¬ 
tion  through  non-planar  layered  media  is  of  practical  interest 

in  the  design  and  construction  of  wideband,  lensed  radomes. 

1-4 

Work  in  this  area  represents  an  exercise  in  the  trans¬ 

mission  of  electromagnetic  radiation  through  a  layered  medium, 
and  most  such  analyses  use  the  approximation  that  the  radia¬ 
tion  penetrates  a  surface  that  can  be  considered  as  "locally 
5-7 

plane."  This  permits  the  use  of  the  "characteristic 

8 

matrices”  for  the  stratified  medium,  and  the  plane  sheet 
transmission  coefficients.  Aside  from  the  manner  of  the 
approximation  is  the  form,  which  is  especially  suitable  for 
calculations  on  a  computer.  It  was  the  purpose  of  this  investi¬ 
gation  to  address  the  question  of  the  accuracy  of  the  "locally 

4 

plane"  approximation,  which  has  been  widely  used  in  the  past,  ' 

7  9-11 

'  and  to  seek  other  approximations  which  might  be  used  to 
more  accurately  determine  the  fields  penetrating  a  non-planar 
layered  medium. 

Most  discussions  of  diffraction  by  a  curved  dielectric 
12-19 

surface  deal  with  its  scattering  characteristics  (e.g., 

backscattering  cross  section)  rather  than  transmission.  Trans¬ 
mission  investigations  through  multilayer  structures  generally 

20 

involve  planes  of  dielectrics  and  these  still  generally 

21  8 
draw  on  the  work  of  Albeles  as  discussed  in  Born  and  Wolf, 

Analysis  of  transmission  through  curved  "electromagnetic 


1  . 

windows"  generally  fall  into  three  categories:  ray-tracing 
ray-tracing  augmented  by  the  plane-sheet  transmission  co- 
9  10  23 

efficients;  '  '  and  "Wave  Spectrum"  surface  integration 

techniques,  including  the  plane  wave  spectrum  (PWS) ^ 

24 

and  the  spherical  wave  spectrum. 

Based  upon  the  results  of  a  library  search  of  previous 

work,  it  was  decided  to  first  investigate  transmission  through 
non-planar,  layered  structures  for  which  exact  solutions  were 
possible,  viz. ,  concentric  cylinders  and  concentric  spheres. 
These  solutions  would  then  lay  a  foundation  for  a  second  stage 
of  investigation  which  would  consist  of  two  parts:  a)  com¬ 
parison  with  results  obtained  from  calculations  using  the 
"locally  plane"  approximation,  and  b)  expansion  of  the  exact 
solutions  in  the  limit  of  large  radius,  hopefully  obtaining 
results  which  would  suggest  suitable  correction  terms  which 
would  serve  to  improve  the  accuracy  of  the  calculations  for 
curved  surfaces. 

Stage  1  work  was  readily  accomplished  by  expanding  on  the 

25 

work  of  Kerker,  et.  al.  Work  in  stage  2,  however  has  not 
been  completed  as  yet.  Attempts  to  expand  the  exact  solutions 
in  the  limit  of  large  radius  for  both  the  cylinder  and*  the 
sphere  have  thus  far  not  yielded  formulations  which  can  be 
readily  implemented  in  transmission  calculations  (i.e.,  appro¬ 
priate  correction  "factors"  or  "terms"  which  could  be  combined 


4 


wish  or  replace  the  "characteristic  matrix"  usually  used  in 
such  calculations) .  Stage  2,  part  (a)  calculations  have  been 
delayed  because  of  a  lack  of  computer  subprograms  to  calculate 
values  for  the  half-integer  Bessel  Functions  and  the  so-called 
"Riccati-Bessel  Functions”.  The  latter  programs  have  been 
written  (in  both  the  BASIC  and  the  FORTRAN  IV  language)  and 
are  included  in  the  appendix  of  this  document.  Calculations 
using  these  programs  are  currently  underway,  but  unfortunately 
have  not  been  completed  for  inclusion  here.  It  is  expected 
that  these  results  will  be  given  in  a  future  communication  and 
publication. 

A  summary  of  the  theoretical  investigations  that  have  been 
carried  out  thus  far  is  given  below. 


II.  Fields  Within  Concentric  Circular  Cylinders 
25  13 

Kerker  and  Matijevic  have  dealt  with  the  case  of 
scattering  by  coaxial  dielectric  cylinders.  Whereas  they  are 
concerned  mainly  with  fields  external  to  the  cylinders  (scatter- 
ing) ,  we  will  consider  the  fields  within.  Consider  plane  elec¬ 
tromagnetic  radiation  traveling  in  x-direction,  incident  upon 
two  coaxial  cylinders.  Let  a*radius  of  the  inner  cylinder 
which  has  index  of  refraction  irij?  let  b*radius  of  outer  cylinder 
of  index  letmQ=index  of  the  space  surrounding  the  cylinders. 


Incident 

Radiation 


AH  indices  of  refraction  will  be  taken  as  real  (non-lossy 
materials) .  The  field  vectors  can  be  related  to  solutions 
of  the  scalar  wave  equation  in  circular  cylindrical  coordi¬ 


nates. 


v*  u  +  u  = 


where  we  have  suppressed  the  harmonic  time  dependence  factor 
This  becomes 

J-  A-(r  -L  ^  (2) 

r  IrV  r*  Ti* 


which 


M  =  R (r)  ©(9)  Z(l). 


This  yields  three  ordinary  differential  equations: 

‘w.1  © o 

(5) 


Where  h  kQsin  ,  with  0  equal  to  the  angle  between  direction 
of  propagation  of  the  incident  radiation  and  a  normal  to  the 
z-axis.  Since  we  are  initially  assuming  the  incident  radiation 
is  along  the  x-axis,  0=0  (0*0  for  oblique  incidence).  For  the 
case  where  the  E-vector  is  incident  parallel  to  the  z-axis. 


the  potentials  in  the  three  regions  are  given  by: 

r>b :  if*,  f 


27 


(7a) 


b>r>a- 


(7b) 


7 


r<cL 


u 


(x) 


a  r  K  Tx(^tkr) 

y\z  -*o  l 


For  the  H-vector  parallel  to  the  z-axis,  we  have 

r>  b :  X J>6*okr)“'d^  Hoofer)  ? 

V®.  |FjA^J^>r)-^W^,fer)J 

Vk  -  -  u>  V. 


(8a) 


(8b) 


r<c l  : 


where 


;f  J.MOj ' 


(8c) 


F  =(-l)nein9 
n 

J  (x)=Bessel  Functions  of  the  first  kind 
n 

Y  (x)=Bessel  Function  of  the  second  kind, 
n 

H  (x)=  J  (x)  -  iY  (x)=Hankel  function  of  the  second  kind 
n  n  n 


and  the  A 


1°)  a  (1)  A_  (1)  A_  (2)#  B__  (D  #  B_  (2)  are 


n  ;  n  -  '  n  "  n  n  -  n 

coefficients  whose  values  are  to  be  determined  from  the 

boundary  Conditions.  Applying  the  boundary  conditions  at  the 

interfaces,  viz. .  Vv\Uj  Vv\  all  continuous 

at  the  boundaries,  we  get  eight  equations: 


8 


vh.  ^5>v,  +■  3~»/S*.)  <9a) 


<0 


0) 


^  >*0  H>v 6*«*i)  "  C>V\^  J*(VH#rf,)  (9b) 


f‘)  7  T  1 


-  ^  V,  H,/v$J)  +  B^'vvm  J*/*,*^  -  3/*,  y  J  =  O 

(9c) 

-  C  v  hM*  tk  j:k^  o 


(96) 


&-[?  v*\  h>\(vk,<)  -  cf^  *?  +  A*  vv»,1  Jv/n*,)  =  "Vo  ^v*«^  ^ 

4-1°^,  M-i  +  A^'V,  3^6v<!)  =  ^o-^vi  (^o^i) 

( 9  f ) 

H*6v04  A^  Vvfj^fvw.^-A^Wli 


(9g) 


-  $■*  "W,  )4>\  (vH^i)+  Av\  ^V>  (w,rft)  ^  6^"*  ^l)zO 


(D 


<  *  k  l> 

<*2.  *  feo- 


(9h) 


The  incident  field  can  be  expanded  in  terms  of  the  Jn  to 


yield:  _  c 


£6f  JJWc 


C  H  9 


(10) 


“M  =  -  *0 


where  kQ=  wave  number  in  the  region  outside  the  cylinders. 


-  9  - 


V 


we  wili.  ne  interested  m  tne  nelds  that  penetrate  into 

re  cylinder,  we  must  evaluate  the  A  3  Usincr 

n  n  3 

(9e-9h),  we  can  write: 

VoHv»(yKo»0  -V|H>6*,e0'  V,  Vh.T*6*#»0 

o  o 

O  ~  x\.Hvv  6*,°6)  V.t  Ts  fv<)  O 


Vvs*  Hv\6*«*0  “  Wv\(vs,a/,)  V,  T^(v4,^i)  O 
V*.#  Hvs^Vo^t)  —  "Vi  ^vs  6*t*0  V,  T*  6*i^<)  O 

o  K*  T*6*.*i)  -\V*0s*O> 


>*«  Hvs(va0^,) 
Vo  6***0 


V,0  Hv\(vv0s(,) 

Vo  Hm(,**'*s0 


“  V( Hm 6*i*0  w,T*6*,«0  v»TwO*o*ti^ 

-  v**  9m  (v*0  v'  3m  6*,®0  v.  T*  (>*■  o *0 

-  VjMw  (\*,*0  W\,?x(v*,rfO  ° 

—va,  Hm  6*t®0  v,' T*  6*w0  o 
_ _ _  » 

“  V|  WjTvsfWjO^  Q 

-  fi*A  (v,o<|)  Vw,  T*  o 

-  V^  K  a  (v,*0  V,  Tv.  fv,  o^O  -  vk  J*  6»i*0 

-  Va .'  V,' Tvs (v^O  T* 6*v*4 


where  we  have  used  ctx  =k#b  .  1 

with  k_  eaual  to  the  wave  number  of  rhe  radiarion  in  the  real on 
r>b.  Eqs.  7c  and  8c  can  then  be  evaluated  to  yield  the  plane 
polarized  components  of  the  incident  beam.  An  arbitrary 
elliptically  polarized  wave  may  be  formed  by  linear  super¬ 
position.  The  extention  to  oblique  incidence  is  made  by  . 

retaining  the  "h"  term  in  Eqs.  4  and  6. 

• 

III  -  Fields  Within  Concentric  Spheres  » 

2  5 

Using  the  notation  of  Kerker,  the  solution  for  diffrac¬ 
tion  by  concentric  dielectric  spheres  is  given  in  terms  of 

"Hertz  -  Debye  potentials"  'TT,  and  fTj.  ,  which  satisfy 

28 

the  wave  equation  in  spherical  coordinates 


r  \  it-1- 

Let  it  *  R6-)  4?^- 

By  separation  of  variables,  this  yields  the  three  ordinary 

differential  equations: 


( 14) 


(  15) 


-  11  - 


(  16) 


"  Bfc  ■  m  *  m 


_i-  4£/st~e 


^(■A+f)  - 


S^S 


o 


(17) 


dr 


+-r 


(18) 


where  n  is  an  integer,  and  m=-n , . . . , o, . . . , +n . 

The  solutions  of  the  radial  equation  (16)  are  Riccati- 
Bessel  functions: 


‘fcCh-)  .  jfT  T 


(19) 


(20) 


iJhere  the  Jn+^  and  Yn+Jj  are  the  half-integer  Bessel  Functions 
of  the  first  and  second  kind.  Also, 

.(*) 


J  *  <k  +  i  X,  •  IW  ^Skr) 


where 


The  solutions  to  Eq.  17  are  the  associated  Lengendre 


=Hankel  functions  of  the  second  kind. 


polynomials. 


©  --  efi'ks «) 


-  12  - 


* 


(21) 


r  it 


j_  Te"-1  fill  y(k  r)  K'h^s) 

i*x  >{vvm)  j 


'J  M  ~  I 


where  we  have  taken  the  incident  radiation  as  plane-polarized 
along  the  z-axis.  The  waves  in  the  three  regions  are  obtained 
from  the  potentials:  . 


(J\ 


k3  vl'  ^ 


(23b) 


a<  r<  b : 


rr^  -  ^r-I  P»%«? )<W 


(24a) 


*U 

rrr,: 


r«t 


(24b) 


™-v*  -  ^c-*)  c«  * 


"K=l 

aO 


(25a) 


-iif-s-  u-a«v  < 


(25b) 


** 
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We  impose  one  boundary  conditions  an  n he  surfaces  r 
r=o  which  yield  a  sec  of  eight  equations: 

^Gm)  +<JL  +- 

vm  ^6*^)  4  V,fx  X#*  (*»*)  4  *)  =  O 

^<w  ^6*v«*)  +V,1  fxXjW*)  4  Ki  s  0 

fv!  No>)  +-L  f  K  U, 

^  kt/-i) = vwx 


where 


k,  b  =  ***>/*, 

^EJej(x  r?ir<a./7j 

“VHi  S  "Ai/vn  >  =  relative  index  of  refraction 

"  f  *  of  core 

=  *VV*/vv  j  =  relative  index  of  refraction  of 
s  wavelength  in  outer  medium 


=a  and 

(26a) 

(26b) 

(26c) 

( 26d) 

( 26e) 

( 26  f ) 

( 26g) 

( 26h) 

layer 
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Since  we  are  interested  in  the  fields  in  the  core,  we  need 


coefficients  . 

h  .  These  car. 

be  obtained 

from  £ 

| 

O 

c  I 

^  » 

t 

j 

C 

o 

*4<W) 

few 

tk) 

•K(^) 

t'M 

o 

v,  tM 

o 

x*(  *o>) 

o 

Wk.  (»«>}) 

0 

o 

o 

W  x«  (’*«•<') 

0 

O 

r;^ 

(>»■$) 

&(■*) 

v,  xK*) 

0 

*sJvv<) 

0 

ViXjfw.-J) 

o 

o 

fc/») 

(27) 


(28) 
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The  coefriciencs 


25  are  thus  determined  and  the  exact 


m  Eos. 

fields  in  the  core  can  be  evaluated  using  the  spherical 

30 

dir.ate  expressions: 


:  .  ±_  _ _  Ji/rlO 

*  V"  dri$  i 'Sue 


c  _  T(r  0 

'rive  7FJ7- 


_  / 6",0 
*  16~ 


H 


-  cW6-L_  Z£r,f)  + 


i 

r  TFT©- 


_j _  sVjO 

v-£~«  '}r  J  <f~ 


-  17  - 


. 11 

-  v^vV> 


coor- 

(29a) 

(29b) 

(29c) 

(29d) 

(29e) 

( 29f ) 


IV.  Asymptotic  Expansions 

For  both  the  cylindrical  and  spherical  cases  described 
above,  asymptotic  expansions  of  the  functions  have  thus  far 
not  proved  fruitful  (the  algebra  and  manipulation  of  the 
determinants  is  rather  unwieldly) .  The  approximations  used 
were  that  b-a«a  and  a  very  large.  This  permits  the  use  of 
the  asymptotic  forms  of  the  functions  and  the  Taylor  expansion 
of  the  derivatives.  This  approach  has  been  temporarily 

suspended  and  will  be  pursued  at  a  later  time. 

V.  Calculations  Using  the  Exact  Solution  for  Concentric  Spheres 
This  evaluation  was  attempted  first  because  it  was  felt  it 

would  offer  the  severest  test  of  the  "locally  plane"  approximation 

since  it  is  three-dimensional.  Computer  programs  to  evaluate 

the  Riccati-Bessel  functions  (or  the  spherical  Bessel  Functions 

of  the  first  and  second  kinds,  J  ,  Y  ,  to  which  they  are  related 

n  n 

by 

or  the  half-integer  Bessel  Functions,  y^  to  which  the^X* 

are  related  bjr  Eqs.  19  and  20)  were  not  available.  It  was  thus 
necessary  to  write  these  programs,  using  several  recursion 
techniques  described  in  the  literature. ^  ^  Documented  versions 
of  these  programs  written  in  FORTRAN  IV  and  BASIC  are  included 


m  the  appendix.  Calculations  using  these  programs  and  equa¬ 
tions  27-29  are  currently  underway.  These  results  will  be 
compared  with  values  obtained  using  the  "locally  plane" 
approximation.  As  originally  proposed,  similar  calculations 
will  be  done  for  the  concentric  cylinder  case. 


IV. 


CONCLUSIONS 


A  library  search  of  materials  accessible  to  me  indicates 
that,  thus  far,  there  has  been  no  systematic  investigation 
(neither  theoretical  nor  experimental)  of  the  "locally  plane" 
approximation  to  transmission  of  electromagnetic  radiation 
through  non-planar  layered  media.  A  theoretical  investiga¬ 
tion  of  the  transmission  into  layered  cylinders  and  spheres 
was  done,  and  calculations  based  upon  these  results  are 
underway.  Computer  programs  in  FORTRAN^  IV  and  BASIC  have  been 
written  to  generate  Riccati-Bessel  Functions  (and  spherical 
and  half-integer  Bessel-Functions)  for  use  in  these  calcula¬ 
tions.  The  investigation  thus  far  has  given  this  author  some 
indication  of  why  such  investigations  have  apparently  not  been 
done  before. 

It  is  expected  that  this  work  will  continue  through  the 
summer  (a  recent  collaborator  of  the  author  has  expressed 
some  interest  in  joining  the  investigation),  and  further 
results  should  be  available  for  submission  for  publication 
in  the  Fall  1980.  It  is  expected  that  these  results  will  give 
an  indication  of  the  theoretical  accuracy  of  the  locally  plane 
approximation,  and  suggest  an  appropriate  experimental  study. 
Improved  approximations  to  the  locally  plane  approximation  are 
still  the  ultimate  goal. 
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SUBROUTINE  HBES1(X,N,  DELTA, RES,  IERR) 

C* 

C* . .  RBES1 :  SUBPROGRAM  TO  CALCULATE  THE  RICC ATI -BESSEL  FUNCTION  OF  THE  1st  KIND. 
C*  THIS  PROGRAM  USES  THE  SUBPROGRAM  BESSP1. 

C* 

C*.. Reference:  Stegun  and  Abramowitz,  Eq.  10.3*1. 

C*  (See  SUBROUTINE  BESSP1  for  complete  references.  Note  also  the  comments  there 
C*  re  the  recursive  nature  of  the  calculations  and  how  to  save  computing  time.) 

C* 

C*..Dr.  Louis  A.  DeAcetis,  February  1980. 

C« 

C* - 

C* 

C***Forra  of  the  CALL  statement:  ” 

C* 

C«  CALL  RBES1(X,N, DELTA, RES, IERR) 

C» 

C*  Input  values: 

C*  X  =  value  of  the  argument;  range:  0  <=  X  <=  100 

C*  (This  must  be  a  DOUBLE  PRECISION  value.) 

C*  N  =  order;  must  be  an  integer  in  the  range:  0  <=  N  <=  100 

C*  DELTA  =  Desired  accuracy  expressed  as  the  fractional  error. 

C» 

C*  Output  values: 

C*  RES  s  results  (value  of  the  Riccati-Bessel  function) 

C«  (This  is  a  DOUBLE  PRECISION  quantity) 

C*  The  mantissas  of  values  of  RES  with  magnitudes  of 

C*  less  than  1.D-70  have  no  significance. 

C*  IERR  =  error  code:  0  if  no  error 

C#  1  if  X,  N,  or  DELTA  out  of  range  (see  above). 

C*  2  if  accuracy  (DELTA)  not  achieved  or  in  doubt. 

C* - 

DOUBLE  PRECISION  X.RES 

C. . .CASE  FOR  X=0,  - >  RES=0. 

IERR=0 

IF(X.NE.O.DO)  GO  TO  100 
<„RES=O.DO 
RETURN 

C... HANDLE  CASE  OF  N=0  HERE.  USE  BESSP1  FOR  OTHER  CASES: 

100  IF(N.NE.O)  GO  TO  110 
RES=DSIN(X) 

RETURN 

110  CALL  BESSP1(X,N, DELTA, RES, IERR) 

RES=X*RES 

RETURN 

END 

0*00000000000000000000000000000000 


fj  O  CJ  fj  n  n  cj  ^  o  o  CJ  CJ  o  ^  o  n  ^  cj  ^  n  CJ  cj  ^  (J  CJ 


C* 


SUBROUTINE  HBES2(X,N,  RES,  IERR) 

..HBES2:  SUBPROGRAM  TO  CALCULATE- THE  RICCATI-BESSEL  FUNCTION 
OF  THE  2nd  KIND.  THIS  PROGRAM  USES  THE  SUBPROGRAM  BESSP2. 

..Reference:  Stegun  and  Abramowitz,  Eq.  10.3. T. 

(See  SUBROUTINE  BESSP2  for  complete  references.  Note  also  the  comments  there 
re  the  recursive  nature  of  the  calculations  and  how  to  save  computing  time.) 

..Dr.  Louis  A.  DeAcetis,  February  1980. 


••Form  of  the  CALL  statement: 

CALL  RBES2(X , N , RES , IERR) 

Input  values:  » 

X  =  value  of  the  argument;  range:  0  <=  X  <=  100 

(This  must  be  a  DOUBLE  PRECISION  value.) 

N  =  order;  mu3t  be  an  integer  in  the  range:  0  <=  N  <=  100 

Output  values: 

RES  a  results  (value  of  the  Riccati-Bessel  function) 

(This  is  a  DOUBLE  PRECISION  quantity.  Tests  indicate 
results  are  accurate  to  better  than  7  digits,  except 
for  values  >  1.D70  in  magnitude.) 

IERR  =  error  code:  0  if  no  error 

1  if  X  or  N  out  of  range  (see  above). 


DOUBLE  PRECISION  X,RES 

C. . .CASE  FOR  X=0,  — >  RES=-1  if  N=0;  RES=  -infinity  if  NX) 

IERR=0 

IF(X.NE.O.DO)  GO  TO  100 
RES=- 1 . DO 
IF(N.EQ.O)  RETURN 
RES=- 1 . D70 
RETURN 

C... HANDLE  CASE  OF  N=0,  X>0  HERE.  GO  TO  BESSP2  FOR  OTHER  CASES: 

100  IF(N.NE.O)  GO  TO  110 
RES=-DC0S(X) 

RETURN 

110  CALL  BESSP2(X,N, RES, IERR) 

RES=X*RES 

RETURN 

END 

COOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO 
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SUBROUTINE  RBES1P(X,N,DELT1,H£S,IERR) 

C* 

C*..RBES1P:  SUBPROGRAM  TO  CALCULATE  THE  DERIVATIVE  OF  THE  RICCATI-BESSEL 
C*  FUNCTION  OF  THE  1st  KIND. 

C*  THIS  PROGRAM  USES  THE  SUBPROGRAM  BESSP1. 

C • 

C*.. Reference:  Stegui  and  Abramowitz,  Eqs.  10.3.1,  10.1.21. 

C*  (See  SUBROUTINE  BESSP1  for  complete  references.  Note  also  the  comments  there 

C*  re  the  recursive  nature  of  the  calculations  and  how  to  save  computing  time.) 

C» 

C*..Dr.  Louis  A.  DeAcetis,  February  1980. 

C* 

C» - 

c* 

C***Form  of  the  CALL  statement: 

C« 

C*  CALL  RBES1P(X,N, DELTA, RES, IERR) 

C« 

C*  Input  values:  • 

C*  X  =  value  of  the  argument;  range:  0  <=  X  <=  100 

C»  (This  must  be  a  DOUBLE  PRECISION  value.) 

C*  N  =  order;  must  be  an  integer  in  the  range:  0  <=  N  <=  100 

C*  DELTA  =  Desired  accuracy  expressed  as  the  fractional  error. 

C« 

C*  Output  values: 

C*  RES  s  results  (derivative  of  the  Riccati-Bessel  function) 

C»  (This  is  a  DOUBLE  PRECISION  quantity) 

C#  The  mantissas  of  values  of  RES  with  magnitudes  of 

C*  less  than  1.D-70  have  no  significance. 

C*  IERR  =  error  code:  0  if  no  error 

C*  1  if  X,  N,  or  DELTA  out  of  range  (see  above). 

C*  2  if  accuracy  (DELTA)  not  achieved  or  in  doubt. 

C« - 

DOUBLE  PRECISION  X,RES,RN,RNLES1 
C... HANDLE  CASE  OF  N=0  HERE;  USE  BESSP1  FOR  OTHER  CASES: 

IF(N.NE.O)  GO  TO  110 

IERR=0 

RES=DC0S(X) 

RETURN 

110  CALL  BESSPKX.N, DELTA, RN, IERR) 

IF(IERR.NE.O)  PRINT  10,  N,X,IERR 
10  FORMATCO  ERROR  IN  SUBPROGRAM  RBES1P  WITH:’,/, 

:  '  Nr  *  ,14,  * ,  X  r  ♦,  D20.H,’;  IERR  r',  12) 

CALL  BESSPKX.N-1, DELTA, RNLES1, IERR) 

RESrX*RN  -  DFL0AT(N)»RNLES1 

RETURN 

END 

COOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO 


c* 

SUBROUTINE  RBE52P(X  ,N, RES .IERR ) 

C* 

C*..HBES2P:  SUBPR0G1 DM  TO  CALCULATE  THE  DERIVATIVE  OF  THE  RICCATT-BESSEL 
C*  FUKTTON  CF  THE  2nd  KIND. 

C*  THIS  PROGRAM  USES  THE  SUBPROGRAM  BESSP2. 

C* 

C*.. Reference:  Stegm  and  Abranowitz,  Eqs.  10.3.1,  10.1.21. 

C*  (See  SUBROUTINE  BESSP2  for  complete  references.  Note  also  the  comments  there 
C*  re  the  recursive  nature  of  the  calculations  and  how  to  save  computing  time.) 

C« 

C*..Dr.  Louis  A.  DeAcetis,  February  1980. 

C* 


C***Forra  of  the  CALL  statement: 

C*  • 

C«  CALL  RBES2P(X , N , RES , IERR ) 

C* 

C*  Input  values: 

C*  *  X  s  value  of  the  argunent;  range:  0  <  X  <=  100 

C*  (This  must  be  a  DOUBLE  PRECISION  value.) 

C*  Ns  order;  must  be  an  integer  in  the  range:  0  <=  N  <=  100 

C* 

C*  Output  values: 

C*  RES  =  results  (derivative  of  the  Riccati-Bessel  function) 

C»  (This  is  a  DOUBLE  PRECISION  quantity) 

C*  The  mantissas  of  values  of  RES  with  magnitudes  of 

C*  less  than  1.D-70  have  no  significance. 

C*  IERR  s  error  code:  0  if  no  error 

C*  1  if  X  or  N  out  of  range  (see  above). 

DOUBLE  PRECISION  X.RES.RN.RNLESI 
C... HANDLE  CASE  OF  N=0  HERE;  USE  BESSP2  FOR  OTHER  CASES: 

C...TEST  X.  MUST  BE  0  <  X  <=  100. 

IF(X.GT.0.D0. AND.X.LE. 1 .D2)  GO  TO  100 

IERR=1 

RETURN 

100  IFCN.NE.O)  GO  TO  110 
IERR=0 
RES=DSIN(X) 

RETURN 

110  CALL  BESSP2(X , N, RN, IERR) 

IF(IERR.NE.O)  PRINT  10,  N,X,IERR  • 

10  F0RMAT(’0  ERROR  IN  SUBPROGRAM  RBES2P  WITH:’,/, 

:  *  Ns  MR,',  X  s  •,  D20.14,*;  IERR  12) 

CALL  BESSP2(X,N-1 ,RNLES1 ,IERR) 

RESsX»RN  -  DFL0AT(N)*RNLES1 

RETURN 

END 

COOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO O 


<2{3<2‘2‘3£2‘2f2cic2‘2****f2<2c2<2‘J<2<2‘2<at2<2‘2<atici‘2<i  ? 


SUBROUTINE  3ESF1Tr(X,N,EELIlTHES,IERR> 

..BESFR1:  SUBPROGRAM  TO  CALCULATE  THE  FRACTIOMAL-BESSHL  FUNCTION 
OF  THE  1st  KIND.  THIS  PROGRAM  USES  THE  SUBPROGRAM  BESSP1. 

..Reference:  Stegui  and  Abraroowitz ,  Eq.  10.1.1. 

(See  SUBROUTINE  BESSP1  for  complete  references.  Note  also  the  comments  there 
re  the  recursive  nature  of  the  calculations  and  how  to  save  computing  time.) 

..Dr.  Louis  A.  DeAcetis,  February  1980. 


‘•Form  of  the  CALL  statement: 

CALL  BESFR1(X,N, DELTA, RES, IERR) 

Input  values: 

X  =  value  of  the  argument;  range:  0  <=  X  <:  100 

(This  must  be  a  DOUBLE  PRECISION  value.) 

N  =  integer  part  of  the  order.  Yields  the  order:  N  ♦  1/2 
must  be  an  integer  in  the  range:  0  <s  N  <=  100 
DELTA  s  Desired  accuracy  expressed  as  the  fractional  error. 

Output  values: 

RES  =  results  =  value  of  the  Fhactional-Bessel  function,  J<N+.5>(X) 
(This  is  a  DOUBLE  PRECISION  quantity) 

The  mantissas  of  values  of  RES  with  magnitudes  of 
less  than  1.D-70  have  no  significance. 

IERR  s  error  code:  0  if  no  error 

1  if  X,  N,  or  DELTA  out  of  range  (see  above). 

2  if  accuracy  (DELTA)  not  achieved  or  in  doubt. 

DOUBLE  PRECISION  X,RES,SQR2PI 
C...SQR2PI  =  SQRT(2/PI) 

DATA  SQR2PI/. 7978845608028654/ 

C. . .CASE  FOR  X=0,  — >  RES=0. 

IERRsO 

IF(X.NE.O.DO)  GO  TO  100 

RESsO . DO 

RETURN 

C... HANDLE  CASE  OF  N=0  HERE.  USE  BESSP1  FOR  OTHER  CASES: 

100  IF(N.NE.O)  GO  TO  110 

RES  a  SQR2PI*DSIN(X)/DSQRT(X) 

RETURN 

110  CALL  BESSP1(X,N, DELTA, RES, IERR) 

RESsSQR2PI *DSQRT ( X ) *  RES 

RETURN 

END 

COOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO 
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c*  •<* 

SUBROUTINE  BESFBtl ,  N  f  RES.lEHf) 

C* 

C*.  .BE5FR2:  StfflPROCRJKTO  CALCDtflE  THE  FRACTTOttL-BESSEL  FUNCTION 
C*  OF  T!E  2nd  KIND.  THIS  PR OCR**  USES  THE  SUBPROGRAM  BESSP2. 

C* 

C*.. Reference:  Stegurr  and  Abraoowitz,  Eq.  10.1.1. 

C*  (See  SUBROUTINE  BESSP2  for  complete  references.  Note  also  the  comments  there 
C*  re  the  recursive  nature  of  the  calculations  and  how  to  save  computing  time.) 

C* 

C*..Dr.  Louis  A.  DeAeetis,  February  1980'. 

C» 

C« 

C***Form  of  the  CALL  statement: 

C» 

C»  CALL  BESFR2(X ,N, RES, IERR) 

C» 

C*  Input  values: 

C*  X  =  value  of  the  argument;  range:  0  <=  X  <=  100  - 

C*  (This  must  be  a  DOUBLE  PRECISION  value.) 

C*  Ns  integer  part  of  the  order.  Yields  the  order:  N  ♦  1/2; 

C*  must  be  an  integer  in  the  range:  0  <=  N  <s  100 

C* 

C*  Output  values: 

C*  RES  s  results  s  value  of  the  Fractional-Bessel  function,  Y<N«-.5>(X) 

C»  (This  is  a  DOUBLE  PRECISION  quantity) 

C*  The  mantissas  of  values  of  RES  with  magnitudes  of 

C*  less  than  1.D-70  have  no  significance. 

C*  IERR  s  error  code:  0  if  no  error 

C*  1  if  X,  N,  or  DELTA  out  of  range  (see  above). 

C* - 

DOUBLE  PRECISION  X,RES,SQR2PI 
C...SQR2PI  s  SQRT(2/PI) 

DATA  SQR2PI/ .7978845608028654/ 

C. . .CASE  FOR  X=0,  - >  RESs-INFINITY  FOR  EVEN  N,  ♦INFINITY  FOR  ODD  N 

IERRsO 

IF(X.NE.O.DO)  GO  TO  100 
RESs-1 ,D70 

IF((N/2)*2-N  .NE.O)  RESs-RES 

RETURN  .♦ 

C... HANDLE  CASE  OF  N=0  HERE.  GO  TO  BESSP2  FOR  OTHER  CASES: 

100  IF(N.NE.O)  GO  TO  110 

RES  s  -SQR2PI*DC08(X)/DSQRT(X) 

RETURN 

110  CALL  BESSP2(X,N, RES, IERR) 

RESsSQR2PI «DSQRT(X )»RES 

RETURN 

END 

COOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO 


c* 

SOBBCOng  3ESSP1  CX,J,EELT.t,BES,IEBR) 

C* 

C* . .  BESSP1 :  SUBPROGRAM  TO  CATJCIILiTE  THE  SPHERICAL  BESSEL  FUNCTION  OF  THE 

c*  iat  kind  i  (x). 

C*  n 
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C***Form  of  the  CALL  statement: 

C» 

C«  CALL  BESSP1CX.N, DELTA, RES, IERR) 

C* 

C*  Input  values: 

C*  X  =  value  of  the  argument;  range:  0  <=  X  <=  100 

C»  (This  must  be  a  DOUBLE  PRECISION  value.) 

C*  Ns  order;  must  be  an  integer  in  the  range:  0  <=  N  <=  100 

C*  DELTA  s  Desired  accuracy  expressed  *«s  the  fractional  error. 

C* 

C*  Output  values: 

C*  RES  s  results  (value  of  the  spherical  Bessel  function) 

C»  (This  is  a  DOUBLE  PRECISION  quantity) 

C*  The  mantissas  of  values  of  RES  with  magnitudes  of 

C*  less  than  1.D-70  have  no  significance. 

C«  IERR  *  error  code:  0  if  no  error 

C*  1  if  X,  N,  or  DELTA  out  of  range  (see  above). 

C*  2  if  accuracy  (DELTA)  not  achieved  or  in  doubt. 


..Because  of  the  recaratve  nature  of  the  calculations  dooe  by 
this  SUBROUTINE,  sane  calculated  values  with  ID5  are  stored  in  an 
array  containing  all  of  the  values  of  the  spherical  Bessel 
function  jn(X)  up  to  the  order  N  that  was  supplied  in  the 
CALL  statement.  For  a  given  value  of  X,  if  sore  than  one  order 
is  to  be  used,  call  this  SUBROUTINE  first  with  the  largest 
value  of  N  to  be  used.  Subsequent  calls  can  then  extract  values 
from  the  array  without  recalculation,' saving  computing  time. 

A  new  calculation  is  initiated  whenever 

(a)  a  new  value  of  X  or  DELTA  is  supplied,  or 

(b)  N  is  larger  than  the  previously  used  maxiraun  value. 

DOUBLE  PRECISION  VALUES( 102) ,X, RES, XOLD, RESOLD 
DOUBLE  PRECISION  A,B,C,D,E,F,T1 ,T2,DS,DC 
DATA  C,X0LD,N0LD/0. ,0. ,0/ 

IERRsO 
RESsO.DO 

IF(N.LT.O.OR.N.GT. 100.OR.X.LT.0. .0R.X.GT.1.D2)IERRs1 
IFCDELTA.LT. 1.E-13)  IERR=1 
IF(IERR.NE.O)  RETURN 
C... HANDLE  CASE  WHERE  X=0  AND  N>0:  RESsO 
RESsO.DO 

IF(X.EQ.O. .AND.N.NE.O)  RETURN 
C...X=0,  N=0:  RES=1;  X>0,  N=0:  RES=SIN(X)/X 
RESs 1 . DO 

IF(X.EQ.O. .AND.N.EQ.O)  RETURN 
IF(N.GT.O)  GO  TO  100 
RES=DSIN(X)/X 
RETURN 

C... FOLLOWING  FOR  X>0,  N>0  CASE: 

100  IF((X.EQ.X0LD.AND.N.LT.N0LD.AND.N0LD.GT.4) .AND. 

:  (C.NE.1.D0.0R.N.GT.4))  GO  TO  400 
DS=DSIN(X) 

DC=DC0S(X) 

IF(N.GT.4)  GOTO  200 
GO  TO  (110, 120, 130, 140), N 
C...N=4: 

140  RESs(((105.D0/(X*X)-45.D0)/X+X)*DS-(105.D0/(X*X)-10.D0)*DC)/X/X 
RETURN 
C...N*3: 

130  RESs((15.D0/X/X-6.D0)*DS-(15.D0/X-X)#DC)/(X#X) 

RETURN 

# .u-2j 

*120  RESs ( ( 3 . DO/X-X) *DS-3 . D0-DC) /(X*X) 

RETURN 
C. ..Nst: 

110  RESs ( DS/X-DC) /X  v 

RETURN  -  %  * 
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G... USING  PREVIOUSLY  GENERATED  71LUES  WHEN  X  UNCHANGED  AND  KLAST  N 
400  RESsC*VALUES(N+1) 

RETURN 

tOI  RES=VALUES(M-»1 ) 

RETURN 
C...X>0,  N>4: 

200  JX=4 

As(15.DO/<X*X)-6.DO)/(X*X) 

B=((  105.D0/(X#X)-45.D0)/X+X)/(X*X) 

Es(105.D0/(X»X)-10.D0>/(X»X) 

F=(X-15.D0/X)/(X»X) 

160  C=(2»JK*1)#B/X  -  A 
JK=JK+1 

D=(1-2*JK)*E/X  -  F 
T1=C*DS 

T2=D*DC  * 

IF(2*( JK/2)-JK.EQ.O)  T2=-T2 

C...IF  THIS  METHOD  BECOMES  TOO  INACCURATE,  00  TO  ALTERNATE. 

IF(T1#T2.LT.0..AND.ABS(TUT2).LTJ.D-14/DELTA»ABS(T1))  GO  TO  500 

VALUES(JKt-1)sTUT2 

IF( JK.EQ.N)  GO  TO  214 

A=B 

B=C 

F=E 

E=D 

GO  TO  160 

C... ALTERNATE  METHOD—  "MILLER'S  METHOD" 

C  FIND  INITIAL  VALUE  FOR  JM  WHICH  IS  LARGER  THAN  N  OR  X  BUT  <=  102 
500  A=DS/X 

RESOLDsO.DO 

JM=MAX(N*1,DINT(X*.5)) 

170  JMsMIN( JM, 101) 

JMsJM+1 

VALUES(JM)=0.D0 
VALUES( JM-1 )=1 . DO 
DO  555  K=3,JM 
KKrJM+1-K 

555  VALUES(KX)*(2 . D0*DFL0AT(KK-1  )«-3 . TO) *VALUES(KK+1 )/X-VALUES(KK*2) 
CsA/VALUES(1> 

RES=C*VAUJE^(N+1) 

IF( ( DABS( ( RES-RESOLD) / RES) . LT. DELTA ). OR .( DABS( RES) . LT . ID-70 ) ) 

:  GO  TO  215 
RESOLD = RES 

IF( JM.GT. 101 )  GO  TO  220 

JMsJM+4 

GO  TO  170 

C... ACCURACY  IN  DOUBT  OR  NOT  ACHIEVED— 

220  IERR=2 
RETURN 

214  RESaVALUES(N*1) 

.  Cal. DO  x 
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...Form  of  the  CALL  statement: 

CALL  BESSP2(X,N,RES, IERR) 

Input  values: 

X  =  value  of  the  argument;  range:  .  0  <=  X  <=  100 

(This  must  be  a  DOUBLE  PRECISION  value.) 

N  =  order;  must  be  an  integer  in  the  range:  0  <=  N  <=  100 

CXitput  values: 

RES  s  results  (value  of  the  spherical  Bessel  function) 

(This  is  a  DOUBLE  PRECISION  quantity.  Tests  indicate 
that  the  results  are  good  to  better  than  seven  significant 
digits,  except  for  magnitudes  >  1.D70.) 

IERR  *  error  code:  0  if  no  error 

1  if  X  or  N  out  of  range  (see  above). 


C***Because  of  the  recursive  nature  of  the  calculations,  all  orders 
of  the  Aaiction  will  be  calculated  up  to  order  N.  If  more  than 
one  order  is  to  be  used  with  the  same  value  of  X,  then  an  initial 
CALL  to  the  SUBROUTINE  with  the  highest  value  of  N  will  cause  all 
lower  orders  to  be  calculated,  and  these  will  be  used  as  long  as 
X  is  unchanged  and  N  <  highest  value.  This  will  save  some  computing  time. 


D()UBLE  PRECISION  VALUES(102) ,X,RES,XOLD 
DOUBLE  PRECISION  DS.DC 
DATA  X0L8 , NOLD/O . , 0/ 

IERR=1 
RES=O.DO 

IFCN.LT.0.0R.N.GT. 100.0R.X.LT.O. . OR. X. GT. 1.D2) RETURN 
*  IERR=0 

C... HANDLE  CASE  WHERE  X=0:  RES=  -INFINITY=-1 .D70 
RES=-1.D70 

IF(X.EQ:0.D0)  RETURN 
DC=DC0S(X) 

IF(N.EQ.O)  GO  TO  105 
C... FOLLOWING  FOR  X>0,  N>0  CASES: 

100  IF(X.EQ.XOLD.AND.N.LT.NOLD)  GO  TO  400 
DS=DSIN(X) 

IF(N.EQ.I)  GO  TO  110 
C.  N=2: 

*120  VALUES(3)=((X-3.DO/X)*DC-3.DO«DS)/(X»X) 

C...N=1: 

110  VALUES(2)=-(DC/X+DS)/X 
C...N=0: 

105  VALUESi 1 )s-DC/X 

IF(N.GT.2)  GO  TO  200 
GO  TO  250 

C... USING  PREVIOUSLY  GENERATED  VALUES  WHEN  X  UNCHANCED  AND  NCHIGHEST  N 
400  RES=VALUES(N*T) 

RETURN 

C... USING  RECURSIVE  RELATIONSHIP 
200  DO  666  K=3,N 

666  VALUES(K*1)=(2.D0»DFL0AT(K)-1.D0)»VALUES<K)/X  -  VALUES(K-I) 

250  RESsVALUES(N+1 ) 

XOLDsX 

NOLDsN 

RETURN 

END 
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13000  REM  SUBPROGRAM  TO  CALCULATE  RfCCATf /SPHERI CAL  BESSEL  FUNCTIONS 
100Q2  REM  OF  THE  1ST  K1  ND , A Mt>  TH€» «  D£<H  NATIVE.. 

10010  REM.. .DR.  LOUIS  A.  DEACETIS 

10011  REM  PHYSICS  DEPARTMENT 

10012  REM  BRONX  COMMUNITY  COLLEGE/CUNY 

10013- REM  BRONX,  NY  10453 

10014  REM 

10015  REM... LAST  UPOATE:  FEBRUARY  1980 

10016  REM 

10018  REM... THIS  WORK  WAS  SUPPORTED  BY  A  GRANT  FROM  THE  USAF  OFFICE  OF 

10019  REM  SCIENTIFIC  RESEARCH,  #790035. 

10020  REM 

10021  REM. . .REFERENCES: 

10022  REM  HANDBOOK  OF  MATHEMATICAL  FUNCTIONS,  ABRAMOWITZ,  M.  ANO  I .A.  STEGUN, 

10023  REM  DOVER  PUBLIC.,  N.Y.1965,  EQS. 10 . 1. 10, 10. 1. 11, 10. 1 . 19, 10. 1. 21, 10. 3. 1 

10024  REM 

10025  REM  STEGUN,  I.  A.,  AND  M.  ABRAMOWITZ,  "GENERATION  OF  BESSEL  FUNCTIONS 

10026  REM  FUNCTIONS  ON  HIGH  SPEED  COMPUTERS,"  MATHEMATICAL  TABLES  AND 

10027  REM  OTHER  AIDS  TO  COMPUTATION,  VOL.  11,  255-257  (1957). 

10029  REM 

10030  REM  JAHNKE,  E.  AND  F.  EMDE,  TABLES  OF  FUNCTIONS,  DOVER  PUBLI CATI ONS, 

10031  REM  NEW  YORK,  1945,  128. 

10032  REM 

10033  REM  KERKER,  M.,  THE  SCATTERING  OF  LIGHT  AND  OTHER  ELECTROMAGNETIC 

10034  REM  RADIATION,  ACADEMIC  PRESS,  NEW  YORK  1969,  PP.  64-67. 

10036  REM  . 

10033  REM  ALL  VARIABLES  AND  FUNCTIONS  USED  BY/ IN  THESE  ROUTINES  BEGIN  WITH  J. 

10040  REM  DO  NOT  USE  VARIABLE  NAMES  BEGINNING  WITH  J  ELSEWHERE. 

10041  REM  THE  NAMES  PI,  P2I,  HALFPI  ARE  ALSO  USED  (*  PI,  2*PI,  PI/2) 

10042  REM  . 


10044  REM  FIRST  CALL  TO  THESE  SUBPROGRAMS  IS  "GOSUB  10000"  WHICH  INITIALIZES 

10046  REM  ALL  FUNCTIONS  AND  DEFINES  CONSTANTS. 

10047  REM 


10043 

10049 

10050 

10051 

10052 

10053 

10054 

10055 

10056 

10057 

10058 

10059 

10060 
10061 


REM 

REM 

REM 

REM 

REM 

REM 

REM 

REM 

REM 

REM 

REM 

REM 

REM 

REM 


TO  GENERATE  VALUES  OF  R I CCAT I -BESSEL  FUNCTIONS,  GOSUB  10100 
TO  GET  VALUES  OF  DERIVATIVES,  GOSUB  10200 
TO  GET  VALUES  OF  SPHERICAL  BESSEL  FUNCTIONS,  GOSUB  10300. 
INPUT  VALUES: 

JX  *  ARGUMENT  (0  <«  JX  < *  100) 

JN  =  ORDER  (0  <*  JN  <»  100)  ,  * 

DELTA  =»  DESIRED  ACCURACY  EXPRESSED  AS  THE  FRACTIONAL  ERROR. 
OUTPUT  VALUES: 

JS  =  VALUE  OF  SPHERICAL  BESSEL  FUNCTION,  JCJNXJX) 

JR  -  VALUE  OF  R I CCAT I  FUNCTION  OR  ITS  DERIVATIVE. 

(MANTISSAS  OF  VALUES  <  IE-70 
JJ  -  ERROR  CODE  »  0  (NO  ERROR),  - 
-  2  (ACCURACY  NOT 


HAVE  NO  SIGNIFICANCE) 

1  (JN,  JX,  DELTA  OUT  OF  RANGE), 
ACHIEVED  OR  IN  DOUBT) 
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10062  REM  BECAUSE  OF  THE  RECURSIVE  NATURE  OF  THE  CALCULATIONS,  FOR  VALUES 

10 06 4  REM  OF  JN  LARGER  THAN  4,.  ALL  ORDERS  OF  THE  FUNCTION  W  I  LL  BE  CALCULATED 

10056  REM  UP  TO  ORDER  JN  .  THUS,  IF  MORE  THAN  ONE  ORDER  WITH  THE  SAME  JX 

10070  REM  IS  TO  3E  USED,  AN  INITIAL  GOSUS  WITH  THE  HIGHEST  VALUE  OF 

10072  REM  OF  JN  TO  BE  USED  WILL  GENERATE  LOWER  ORDERS  AS  *ELL,  AND 

10074  REM  THESE  WILL  BE  USED  AS  LONG  AS  JX  AND  DELTA  ARE  UNCHANGED  AND 

10076  REM  JN<H I GHEST  VALUE  (THIS  CAN  SAVE  CONSIDERABLE  COMPUTING  TIME). 

10080  REM  - 

10081  REM  ROUTINES  FOR  SINE  AND  COSINE  TO  15  DECIMAL  PLACES  ARE  AT  LINES 

10082  REM  10900  -  10990.  INPUT  ANGLE  IS-AA. 

10084  REM  GOSUB  10900  FOR  SIN(AA)  RETURNED  AS  VARIABLE  SA. 

10086  REM  GOSUB  10950  FOR  COS(AA)  .RETURNED  AS  VARIABLE  CA. 

10088  REM  IF  THE  BASIC  BEING  USED  HAS  DOUBLE  PRECISION  FUNCTIONS,  THEN  OMIT 

10089  REM  THESE  LINES,  AND  CHANGE  "SA"  AND  "CA"  TO  SIN  AND  COS  BELOW. 

10090  REM . . 

10091  DIM  JF(101) 

10092  Pl=3. 141592653589793:  P2 I =6 . 23  3185307179586  :  HALFPI =1. 570796326  794896 
10094  DEF  FNA(A)=A-INT(A/P2I )  *  P2  I  : REM  MAKE  "A"  A  VALUE  BETWEEN  0  AND  2  PI 
10096  DEF  FNT(N,X)=(2*N+3)*JF(N+l)/X-JF(N+2) 

10098  JC=0:  JZ  =  0 :  JY  =  0:  RETURN 

10100  REM  FIND  THE  Rl CCAT I -BESSEL  FUNCTION  OF  THE  FIRST  KIND: 

10102  IF  JN < 0  OR  JXCO  OR  DELTAC1.E-13  THEN  JR  =  0:  JJ  =  1:  RETURN 

10106  JJ  =  0 

10108  IF  JX  =  0  THEN  JR  =  0:  RETURN:  REM  JR=0  IF  JX=0 

10110  IF  JNOO  GOTO  10120 

10115  AA= JX :  GOSUB  10900:  JR=SA:  RETURN 

10120  GOSUB  10300 

10150  JR=JX* JS:  RETURN 

10107  REM  HANDLE  JX,JN=>  0  CASES  HERE:  USE  SPHERICAL  BESSEL  ROUTINE  FOR  OTHERS: 

10200  REM  FIND  FIRST  DERIVATIVE  OF  THE  R I CCAT I -BESSE L  FUNCTION  OF  FIRST  KIND: 

10201  IF  JN<  OOR JX  <  0  OR  DELTAC1.E-15  THEN  JR=0:  JJ  =  1:  RETURN 

10202  JJ=»0:  IF  JNOO  GO  TO  10206 

10204  AA=JX :  GOSUB  10950:  JR=CA:  RETURN 

10205  J  J=*0 

10206  GOSUB  10300:  JT=JS:  IF  JJOO'PRINT  'ERROR  FOR  JN,JX=  ';JN,JX 

10208  JN=JN-1 :  G0SU3  10300:  JN=JN*1 

10209  JR=»JX*JT-JN*JS:  RETURN 

10300  REM  FIND  SPHERICAL  BESSEL  FUNCTION  OF  THE  FIRST  KIND: 

10301  IF  JN< 0  OR  JXCO  OR  DELTAC1.E-13  THEN  JS  =  0:  JJ  =  1:  RETURN 

10305  IF  JX=*0  AND  JNO  0  THEN  JS  =  0 :  RETURN:  REM  JS  =  0  FOR  JX  =  0  EXCEPT  IF  ^  =  0 

10311  IF  JX*JY  AND  JN<=JZ  AND  JZ>4  AND  (JCOl  OR  JN>4)  GOTO  10400 

10312  AA-JX:  GOSUB  10900;  IF  JNOO  GOTO  10316 

10314  IF  JX=0  THEN  JS»1  ELSE  JS=SA/JX 

10315  RETURN 

10316  GOSUB  10950:  ON  JN  GOTO  10320,10330,10340,10350 
10318  GOTO  10360:  REM  JN>=5 

10320  JS*(SA/JX  -  CAJ/JX:  RETURN  :  REM  JN-1 
10330  JS*( { 3/JX  -  JX)*SA  -  3*CA) /( JX* JX ) :  RETURN:  REM  JN=2 
10340  JS-(( 15/( JX*JX)-6)*SA-( 15/JX-JX)*CA)/(JX*JX):RETURN:  REM  JN-3 
10350  JS»((( 105/(JX*JX)-45)/JX>JX)*SA-(105/( JX* JX ) -10 ) *CA) / ( JX* JX > 
v  10352  RETURN  : REM  JN-4 


10360  JA={15/CJX*JX)-6  )/(  JX*JX)  :  JB=(  (105/  ( JX*  JX ) -4  5  )  /  JX+ JX)  /  (  JX*JX  > 

10370  JE*C105/(JX*JX)-1Q)/(JX*JX)  :  JF  =  ( JX-15/ JX) / ( JX*JX ) 

10380  JC*C2*JK+1)*JB/JX  -  JA 
10385  JK=JK+1 

10390  JD=*(  1-2*JK)*JE/JX-JF 

10594  T1*JC*SA:  T2=JD*CA:  I F ( 2* I NT( JK/ 2 ) - JK) =0  THEN  T2— T2:  REM  (-1)|(N+1) 

10596  IF  T1*T2  <0  AND  ABS( T1+T2 ) < 1 . E-14 /DELTA*ABS< T1 )  GOTO  10500 

10597  JF(JK)=T1*T2 

10598  IF  JK-JN  THEN  J3=JF(JN):  JC=1:  GOTO-10600 

10599  JA=JB:  JB=JC:  JF=JE:  JE=JD:  GOTO  10380 

10400  REM  USING  PREVIOUSLY  GENERATED  VALUES  WHEN  JX  UNCHANGED  AND  JN<OLD  JN 

10410  JS=JC*JF( JN) :  RETURN 

10500  PRINT  "USING  MILLER'S  METHOD" 

10502  JA=SA/JX  : REM  FIND  LARGER  VALUE  OF  JX,JN  <100: 

•  10  505  JF=0:JM31NT(-JX*(JX>JN)-JN*(JN>JX)+.5):REM  CHOOSE  JM=MAX( JX, JN ) +1<=100 
10508  JM— JM*( JM<101)-100*( 100<JM) :  JM=JM+1 
10510  JF( JM)=0:  JF(JM-1)=1 
10520  FOR  JK= JM-2  TO  0  STEP  -1 
10530  JF( JK)=FNT(JK, JX):  NEXT  JK 
105 5=0  JC=JA/JF(0):  REM  JC=  SCALING  CONSTANT 
10555  JS=JC* JF ( JN ) 

10571  REM  FINISHED  IF  VALUE  <  IE-70  OR  DESIRED  ACCURACY  ACHIEVED 

10572  I  F(A8S((  JS- JF ) / JS ) <DELTA )  OR  (  ABS(  JSXl.  E  —  7  0  )  GOTO  10600 
10574  JF=JS:  IF  JM<100  THEN  JM=JM+5:  GOTO  10508 

10576  JJ=2 
10580  RETURN 

10600  JZ-JN:  JY=JX :  RETURN:  REM  SAVE  OLD  VALUES  FOR  USE  IN  10400 

10900  REM . **SER I ES  FOR  SIN  &  COS** . 

.10905  REM  INPUT  ANGL E=AA .  OUTPUT:  SA=SIN(AA);’  CA=COS(AA) 

10907  REM, .VARIABLE  NAMES  USED:  AA,  CA,  SA,  SN,  XX  CONSTANTS:  PI,  HALFPI 

10908  REM  FUNCTION:  FNA(X) 

10910  REM. . . S I NE( AA, : 

10915  IF  AA=0  THEN  SA*0:  RETURN 

10920  XX=FNA( AA) :  SN=1:  IF  XX>PI  THEN  SN— 1:  XX=XX-PI 
10925  IF  XX>HALFP I  THEN  XX=PI-XX 

10930  SA-1-XX*XX/110*( 1-XX*XX/ 156 *( 1-XX*XX/21 0*( l-XX*XX/272*( l-XX*XX/342 ) ) ) ) 
10935  SA=XX*( l-XX*XX/6*( l-XX*XX/20*( l-XX*XX/42*( 1-XX*XX/72*SA) ) ) ) 

10940  SA=»SN*SA:  RETURN 
10950  REM . . COS  I NE( AA) : 

10955  IF  AA=0  THEN  CA=1:  RETURN 

10960  XX=FNA( AA) :  SN=1:  IF  XX>PI  THEN  XX=XX-PI :  SN— 1 

10965  IF  XX>HALFP I  THEN  XX-PI-XX:  SN  — SN  , 

10970  CA=1-XX*XX/132*(  1-XX*XX/182*( l-XX*XX/240*( l-XX*XX/306*(  l-XX*XX/380) ) ) ) 
10975  CA*l-XX*XX/2*( 1-XX*XX/12*( l-XX*XX/30*( l-XX*XX/56* ( 1-XX*XX/90*CA) ) )) 
10980  CA-SN*CA:  RETURN 

10990  REM . ** END  OF  SIN  &  COS  ROUTINES** . 
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REM  SUBPROGRAM  TO  CALCULATE  R l CCAT l /SPHER I  CAL  BESSEL  FUNCTIONS 
REM  OF  THE  2ND  KINO  AND  THE  I R  DERIVATIVE. 

REM.. .OR.  LOUIS  A  DEACETI  S 

REM  PHYSICS  DEPARTMENT 

REM  BRONX  COMMUNITY  COLLEGE/CUNY 

REM  BRONX,  NY  10453 

REM 

REM.. LAST  UPDATE:  FEBRUARY  1930 
REM 

REM... THIS  WORK  WAS  SUPPORTED  BY  A  GRANT  FROM  THE  USAF  OFFICE  OF 
REM  SCIENTIFIC  RESEARCH,  #790035 

REM . 

REM. .REFERENCES: 

REM  HANDBOOK  OF  MATHEMATICAL  FUNCTIONS,  ABRAMOWITZ,  M.  AND  I.A.  STEGUN* 
REM  DOVER  PUBLICATIONS,  NEW  VORK  1965,  EQS  .  10 . 1. 12, 10.119, 10 . 3.  1, 10.1 . 2  1 
REM 

REM  STEGUN,  I.A.,  AND  M  ABRAMOWITZ,  "GENERATION  OF  BESSEL  FUNCTIONS 
REM  FUNCTIONS  ON  HIGH  SPEED  COMPUTERS,"  MATHEMATICAL  TABLES  AND 
REM  OTHER  AIDS  TO  COMPUTATION,  VOL.  11,  255-257  (1957). 

REM 

REM  JAHNKE,  E.  AND  F.  EMDE,  TABLES  OF  FUNCTIONS,  DOVER  PUBLICATIONS, 

REM  NEW  YORK,  1945,  126-129. 

REM 

REM  KERKER,  M.,  THE  SCATTERING  OF  LIGHT  AND  OTHER  ELECTROMAGNETIC 
REM  RADIATION,  ACADEMIC  PRESS,  NEW  YORK  1969,  PP.  64-67. 

REM . 

REM  ALL  VARIABLES  AND  FUNCTIONS  USED  BY/IN  THESE  ROUTINES  BEGIN  WITH  Y. 
REM  DO  NOT  USE  VARIABLE  NAMES  BEGINNING  WITH  Y  ELSEWHERE. 

REM  THE  NAMES  PI,  P2 I  AND  HALFPI  ARE  ALSO  USED  ( -  PI,  2*PI,  PI/2) 

REM . 

REM  FIRST  CALL  TO  THESE  SUBPROGRAMS  IS  "GOSUB  10000"  WHICH  DEFINES 
REM  ALL  FUNCTIONS  ANQ  CONSTANTS,  AND  INITIALIZES  VARIABLES 
REM 

REM  TO  GENERATE  VALUES  OF  THE  Y-R I CCAT I  -  BESSEL  FUNCTIONS,  GOSUB  10100 
REM  TO  GET  VALUES  OF  THE  DERIVATIVE  OF  THE  Y-RICCATl'S,  GOSUB  10200 
REM  TO  GET  THE  SPHERICAL  BESSEL  FUNCTIONS  OF  THE  2ND  KIND,  GOSUB  10300 
REM  INPUT  VALUES: 

REM  YX  *  ARGUMENT  (0  <=*  YX  <*  100)  (0  <YX  *<*  100  FOR  DERIVATIVE) 

REM  YN  »  ORDER  (0  <«  YN  <■  100) 

REM 

REM  OUTPUT  VALUES: 

REM  YS  -  VALUE  OF  SPHERICAL  BESSEL  FUNCTION  Y<YN>(YX) 

REM  YR  -  VALUE  OF  THE  Rl CCATI -BESSEL.  FUNCT I  ON  OR  ITS  DERIVATIVE 

REM  (ACCURACY  FOR  YR  AND  YS  BETTER  THAN  5  DIGITS) 

REM  JJ  -  ERROR  CODE  *  0  (NO  ERROR),  -  1  (YN,  YX  OUT  OF  RANGE) 

REM . - . 
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1006 2  REX  BECAUSE  OF  THE  RECURSIVE  MATURE  OF  THE  CALCULATIONS,  ALL  ORDERS 
10064  RB*  OF  THE  FUNCTION  WILL  BE  CALCULATED  UP  TO  ORDER  YN  IF  MORE  THAN 

IQ 066  REM  ONE  ORDER  IS  TO  BE  USED  WITH  THE  SAME  VALUE  OF  YX,  THEN  AN  INITIAL 

10072  REM  G0SU8  WITH  THE  HIGHEST  VALUE  OF  YN  WILL  CAUSE  ALL  LOWER 
10074  REM  ORDERS  TO  BE  CALCULATED  AND  THESE  WILL  BE  USED  AS  LONG  AS  YX 

10076  REM  IS  UNCHANGED  AND  YN<HI GHEST  VALUE  (THIS  CAN  SAVE  CONSIDERABLE 

10078  REM  COMPUTING  TIME). 

10080  REM  . 

10081  REM  ROUTINES  FOR  SINE  AND  COSINE  TO  15  DECIMAL  PLACES  ARE  AT  LINES 

10082  REM  10900  -  10990.  INPUT  ANGLE  IS  AA. 

10084  REM  GOSUB  10900  FOR  SIN(AA)  RETURNED  AS  VARIABLE  SA. 

10086  REM  GOSUB  10950  FOR  COS(AA)  RETURNED  AS  VAR, ABLE  CA. 

10088  REM  IF  THE  BASIC  BEING  USED  HAS  DOUBLE  PRECISION  FUNCTIONS,  THEN  OMIT 

10089  REM  THESE  LINES,  AND  CHANGE  "SA"  AND  "CA"  TO  SIN  AND  COS  BELOW. 

10090  REM . . . 

10091  DIM  YF(101) 

10092-  PL«3. 14159  2653589  783:  P2L=6. 2  8  318  5307179  586 :  HALF P I »1. 570796 326794886 - 

10094  DEF  FNA(  A)  =*A-  I  NT{  A/P2  I  )  *  P2  I  : REM  MAKE  "A"  A  VALUE  BETWEEN  0  AND  2  PI 
10096  DEF  FNT(N,X)=(2*N-l)*YF(N-l)/X-YF(N-2) 

10098  YZ=0 :  YY»0 :  RETURN 

10100  REM  FIND  Rl CCATT I -BESSEL  FUNCTION  OF  THE  SECOND  KIND: 

10103  IF  YN< 0  OR  YX<0  THEN  YR=0:  JJ  =  1:  RETURN 
10105  JJ-0 

10108  IF  YX=OANDYNOOTHEN  YR«-1 . E7Q: RETURN: REM  YR=-  INFINITY  FOR  YX=0  AND  YNOO 
10110  IF  YX=0  AND  YN=0  THEN  YR=-1:  RETURN:  REM  YR=-1  FOR  YX=0  AND  YN=0 
10120  GOSUB  10311:  YR=YX*YS:  RETURN 

10200  REM  FIND  FIRST  DERIVATIVE  OF  R I CCAT I -BE SSEL  FUNCTION  OF  THE  SECOND  KIND: 

10210  IF  YX>0  AND  YX<=100  GOTO  10220 

10212  JJ*1:  RETURN 

10220  JJ  =  0:  IF  YNO  GOTO  10230 

10222  AA=YX:  GOSUB  10900:  YR=SA:  RETURN:  REM  YN=0  »=>  DER I V=S I N( YX ) 

10230  GOSUB  10300:YT»YS:  I F  JJOO  THEN  PRINT’ERROR  IN  DERI  V:  YN,YX=  ';YN,YX 
10232  YN*YN-1:  GOSUB  10300:  YN=YN*1 
10234  YR»YX*YT-YN*YS:  RETURN 

10300  REM  FIND  SPHERICAL  BESSEL  FUNCTION  OF  THE  SECOND  KIND: 

10302  IF  YN<0  OR  YX<0  THEN  YS»0:  JJ=1:  RETURN 

10305  JJ*0:  IF  YX*0  THEN  YS  —  1.E70:  RETURN:  REM  YS—INFINITY  FOR  YX«=0 

10311  l«  YX=YY  AND  YN<-YZ  GOTO  10400 

10312  AA-YX:  GOSUB  10950:  IF  YN*0  GOTO  10345 
10316  GOSUB  10900:  ON  YN  GOTO  10340,10330 

10330  YF( 2) »( ( YX  -  3/YX)*CA  -  3*SA)/CYX*YX):  REM  YN-2 

10340  YF(D— (CA/YX  ♦  SA)/YX:  REM  YN*1 

10345  YF(0)  —  CA/YX:  REM  YN-0 

10350  IF  YN>2  GOTO  10500 

10360  GOTO  10540 
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10400  REM  USING  PREVIOUSLY  GENERATED  VALUES  .WEN  YX  UNCHANGED  AND  YN<OLD  YN 
10410  YS*YF(YN) :  RETURN 

10500  PRINT  ’USING  RECURSION  RELATIONSHIP' 

10510  FOR  YK=3  TO  YN 
10520  YF(YK)=FNT(YK,YX) 

10530  NEXT  YK 
10540  YS=YF( YN) 

10600  YZ=YN :  YY=YX :  RETURN:  REM  SAVE  OLft  VALUES  FOR  USE  IN  10400 

10900  REM . **SERI ES  FOR  SIN  &  COS** . 

10905  REM  INPUT  ANGLE=AA.  OUTPUT:  SA=SIN(AA);  CA=COS(AA) 

10907  REM.. .VARIABLE  NAMES  USED:  AA,  CA,  SA,  “SN,  XX.  CONSTANTS:  PI,  HALFPI 

10908  REM  FUNCTION:  FNA(X) 

10910  REM.. .SINE(AA): 

10915  IF  AA=0  THEN  SA=Q:  RETURN  , 

10920  XX=FNA(AA):  SN-1:  IF  XX>PI  THEN  SN*-1:  XX=XX-PI 
10925  IF  XX>HALFPI  THEN  XX*PI-XX 

10930  SA=1-XX*XX/110*( 1-XX*XX/156*( 1-XX*XX/£10*( l-XX*XX/2 72* C l-XX*XX/34  2)  )  )  ) 
10935  SA=XX*(l-XX*XX/6*(l-XX*XX/20*(l-XX*XX/42*( 1-XX*XX/72*SA) ) )) 

10940  SA3SN*$A:  RETURN 
10950  REM. ..COSINE(AA): 

10955  IF  AA=0  THEN  CA=1:  RETURN 

10960  XX=FNA(AA) :  SN=1:  IF  XX>PI  THEN  XX=XX-P I:  SN=-1 
10965  IF  XX>HALFP  I  THEN  XX=PI-XX:  SN  —  SN 

109  70  CAal-XX*XX/132*(l-XX*XX/182*(l-XX*XX/240*( l-XX*XX/305*( l-XX*XX/380  ) ) )  ) 
10975  CA=l-XX*XX/2*( 1-XX*XX/12*( l-XX*XX/30*( l-XX*XX/56* ( l-XX*XX/9 0*CA) ) )) 
10980  CA=SN*CA:  RETURN 

10990  REM . **END  OF  SIN  4  COS  ROUTINES** . 


